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ADDRESSING  TENSION  INSTABILITY  IN  SPH  METHODS 


1.  INTRODUCTION 

Smoothed  particle  hydrodynamics  (SPH)  is  a  Lagrangian  method  that  requires  no 
spatial  mesh.  SPD  is  really  an  interpolation  method  in  which  particles  can  be  em¬ 
ployed  as  part  of  the  approximation  [1].  (Note  that  in  this  paper,  the  words  element 
and  particle  will  be  used  interchangeably  when  referring  to  SPH  methods).  The  calcu¬ 
lation  of  interactions  among  the  particles  is  based  upon  their  separation  alone.  This  as¬ 
pect  along  with  the  absence  of  a  grid  allows  very  large  deformations  to  be  computed 
in  a  straightforward  fashion.  In  addition,  arbitrary  fracture  surfaces  can  be  opened 
without  requiring  special  considerations  such  as  an  a  priori  knowledge  of  the  fracture 
location  as  in  finite  element  methods  (FEM).  SPH  is  thus  an  appealing  and  valuable 
computational  tool,  especially  for  high  deformation  events  such  as  impact. 

However,  SPH  is  a  maturing  method  and  still  has  a  few  technical  barriers  to  over¬ 
come  before  becoming  a  widely  used  tool  in  computational  mechanics.  Its  biggest 
drawback  is  the  well  know  instability  in  tension  -  premature  fragmentation  of  the  SPH 
grid  in  tension.  Swegle  et  al.  [2]  have  done  a  formal  stability  analysis  that  clearly  dis¬ 
cusses  the  roots  of  tension  instability.  Hicks  et  al.  [3]  have  applied  a  conservative 
smoothing  approach  with  some  success,  but  tension  instability  remains  a  serious  ob¬ 
stacle  in  general.  Other  difficulties  with  SPH  methods  include  the  accurate  calculation 
of  finite  strains,  inhomogeneous  media,  stress  oscillations  or  “sawtooth”  behavior  (see 
[4]),  and  the  application  of  natural  boundary  conditions. 

In  this  paper,  we  present  a  new  unconventional  (UC)  approach  to  address  tension 
instability  (and  also  oscillatory  stresses)  in  which  the  stresses  are  calculated  at  points 
other  than  the  SPH  nodes.  In  addition,  stress  rate  calculations  are  briefly  addressed  and 
specialized  to  one  dimension  (ID)  applications. 
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2.  THE  SPH  METHOD 


The  SPH  method  is  non-local  in  nature  and  based  upon  interpolation  -  see  [1],  [4- 
6].  Consider  the  function  /(x)  and  its  kernel  average  (/(x)  )  expressed  as: 

(/■(?)>  =  (1) 
where  V  is  the  volume,  x  is  the  current  position  vector,  and  W  is  the  kernel  with  a 
width  measured  by  the  parameter  h.  Various  possibilities  exist  for  the  choice  of 
1V(x  —  x',h)  with  the  most  popular  being  the  cubic-b-spline  [4].  The  kernel  W 
should  satisfy  the  following  condition: 

W(x-x',h)  =  5(x-x’)  (2) 

where  5  is  Dirac’s  delta  function.  Normalization  of  W  and  compact  support  (W  is 
zero  outside  a  limited  domain)  are  expressed  as: 

(3a) 

and 

W(x,h)  =  0  (\x\>2h)  (3b) 

Associating  with  a  particle  j,  a  volume 

dV>  =  K  (4) 

P' 

introduces  the  concept  of  particle  mass,  ,  and  density,  .  The  integral  in  (1)  can 

then  be  replaced  by  the  approximation: 

N  j 

(/'(x)>~  ^  (5) 

7=1  P' 

where  N  is  the  total  number  of  particles. 

Dropping  the  (  )  sign,  gradients  of  /(x)  can  be  expressed  as  [4]: 


2 


(6) 


N 


V(/(J))  =-Y  ITL  /(*/) 

J=lp’ 


cW  J 

M  ‘ 


in  which  the  subscripts  denote  component  and  is  the  /  component  of  the  unit  vector 
at  j.  Another  form  for  the  gradient  of  /(x)  that  is  widely  used  but  which  introduces 
additional  approximation  is  [4]: 


N 

V  (fix) )  =  -p  (X)  ^ 

;=i 


/(?) 

.P^  (?) 


(?•') . 


^  ei 


(7) 


3.  THE  SPH  EQUATIONS  (BRIEFLY) 

In  subscript  notation,  the  conservation  of  linear  momentun  is  expressed  as: 


dV 


m 


1  3a 

1  mn 


dt 


(8) 


in  which  is  the  velocity  and  is  the  Cauchy  stress  tensor,  and  the  sununation 
convention  is  implied  by  the  repeated  index  n .  Applying  (6),  (8)  becomes  at  particle  i: 


dt 


P'  J.l  P'  P^’  K 


(9) 


This  form  does  not  conserve  momentum  since  the  force  on  particle  i  due  to  j  is  not  the 
same  as  j  due  to  i.  A  widely  used  alternative  to  (9)  that  conserves  momentum  but  in¬ 
troduces  further  approximation  is  obtained  by  substituting  (7)  into  (8)  to  produce  (no 
sununation  on  i): 


dv[ 


N 


m 


dt 


7  =  1 


G^ 


mn 


G‘ 


-I- 


mn 


(P^)  (pO 


3jc{, 


(10) 
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Artifical  viscosity  is  usually  included  in  the  linear  momentun  as  an  artifical  vis¬ 
cous  pressure.  Benz[l],  Libersky  and  Petschek  [2]  and  Libersky  etal.  [7]  give  details 
for  including  it  into  (10). 

Mass  conservation  is  of  the  form; 


^ 

dt 


(11) 


Introducing  (7)  into  (11)  produces: 

=  -o‘  V^vi  ^ 
*  ,  O'- 

A  widely  used  alternative  to  (12)  is  [4]: 


dt 


The  density  can  also  be  directly  determined  from  (5) 


dW^j 


and  is  of  the  form: 


(12) 


(13) 


N 

p'  =  (14) 

7=1 

Equation  (14)  will  be  employed  in  this  paper. 

In  addition  to  momentum  and  mass  conservation,  energy  conservation  can  also  be 
included  in  the  governing  equations  (see  [1]  and  [2]  for  instance).  In  this  work,  it  has 
not  been  introduced  since  a  linear  elastic  material  will  be  assumed. 
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4.  CONSTITUTIVE  EQUATIONS  AND  STRESS  RATES 
The  current  position  X  (t)  of  the  material  point  X  can  be  expressed  as; 


X  (t)  =  x  +  U(t)  (15) 

where  U  (t)  is  the  current  displacement  as  measured  from  x.  In  index  notation,  the 
rate  of  deformation  (t)  is  of  the  form: 

1  fdV^  dV\ 

Note  in  (16),  that  d^^  is  an  Eulerian  variable  and  the  derivatives  are  with  respect  to 
x^  (not  X^.  The  rate  of  deformation  can  be  put  in  SPH  form  by  directly  applying  ei¬ 
ther  equation  (6)  or  (7).  In  [2],  further  manipulations  and  approximations  are  intro¬ 
duced  to  produce  the  correct  trace  and  to  express  d  in  terms  of  velocity  differences. 
For  ID  applications,  the  use  of  (6)  or  (7)  will  suffice. 

Next  consider  the  velocity  gradient  which  is  defined  by: 

dV 

(17) 

Using  (16),  can  be  divided  into  symmetric  and  skew-symmetric  parts: 

L  =  d 
mn  mn  mn 

where  is  the  rotation  and  given  by 
mn  °  •’ 

™  2  [dx„  dxj 

and  d^^  is  determined  from  (16). 

Having  defined  the  rate  of  deformation,  we  next  seek  an  objective  stress  rate.  The  Jau- 
man  stress  rate  6^^  is  widely  used  [8],  and  is  expressed  as: 


(18) 

(19) 
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(20) 


=  a  -6, .a  ,-0,  o  , 

mn  mn  Ij  ml  Im  nl 


where  and  6^^  are  components  of  the  Cauchy  stress  and  the  rate  of  Cauchy 
stress  (which  is  not  objective).  For  ID  applications,  0^^  =  0  and 
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(21) 


The  Jauman  stress  rate  can  be  related  to  the  rate  of  deformation  via  the  constitutive 
tensor  C.  In  component  form  this  relation  is: 

a  =  C  t  di 

mn  mnlq  Iq 

For  ID  applications  and  a  linear  material,  (22)  reduces  to: 

^11  =  ^1111^11 


where  Cj  j  j  j  =  E  (the  modulus  of  elasticity  and  assumed  to  be  constant),  and 


‘‘ii  ~ 


(24) 


Explicit  integration  of  Oj  j  to  obtain  Gj  j  (0  is  in  the  form: 

<^11  (0  =  +a^^{t-dt)dt  (25) 

where  dt  is  the  current  time  step  and  from  (23) 

G.^.^(t-dt)  =  E  •  d^^{t  -  dt)  dt  (26) 

Because  the  displacement  t/ j  is  the  time  derivative  of  and  E  is  constant,  (26) 
can  be  directly  integrated  in  time  to  yield: 

Gjj  =  Ee^^  (27) 

where  through  (24) 
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(28) 


Equations  (6)  or  (7)  are  then  introduced  into  (28)  to  express  £  j  |  (0  in  a  SPH  form. 
We  note  that  e^j,  which  is  determined  by  integrating  in  time,  is  not  the  well 
known  linear  strain.  Rather  it  is  an  Eulerian  variable  because  the  derivative  in  (28)  is 
with  respect  to  (t)  not  Jcj . 

5.  ADDRESSING  TENSION  INSTABILITY 
Tension  instability  is  a  problem  that  has  longed  troubled  conventional  SPH  meth¬ 
ods  and  greatly  limited  its  application.  In  [2],  a  stabiltiy  analysis  for  ID  is  presented. 
That  analysis  shows  SPH  calculations  to  be  unstable  when: 

a„>0  (29) 

This  means  that  the  standard  SPH  method,  in  which  the  node  and  the  stress  point  are 
both  located  at  the  centroid  of  the  element  (particle),  will  be  unstable  in  tension.  Swe- 
gle  et  al.  in  [2]  also  demonstrate  that  this  tension  instability  can  not  be  corrected  by 
artifical  viscosity. 

The  approach  in  this  work  to  overcome  tension  instability  is  indicated  in  Fig.  1,  in 
which  the  stress  points  (denoted  by  x)  are  computed  at  points  away  from  the  SPH 
nodes  (denoted  by  o  and  located  at  the  center)  as  measured  by  the  distance  R.  The  lim¬ 
its  of  R  are: 

0<R<0.5  (30) 

R  =  0.5  corresponds  to  the  conventional  form  in  which  both  stress  points  are  located 
at  the  centroid  of  the  element,  while  R  =  0  places  the  two  stress  points  at  the  left  and 
right  edges  of  the  element.  In  conventional  SPH  because  the  stress  points  coincide 
with  the  SPH  node,  particle  i  does  not  enter  into  the  calculation  of  G  ^  j ,  since 


<.dx^  > 
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— :  =  0  for  a  symmetric  kernel  when  i  =  j.  In  this  new  approach,  which  we  re- 
^?to  as  unconventional  or  UC,  for  0  <  /?  <  0.5  node  i  will  be  included  in  the  stress 
calculations  for  SPH  element  i.  Thus  stress  calculations  from  the  unconventional  ap¬ 
proach  should  be  more  accurate  since  node  i  will  be  included  at  the  stress  points 
and  (see  Figure  1).  Also  the  UC  method  should  help  to  lessen  or  eliminate  oscilla¬ 
tion  or  the  “sawtooth”  effect  evident  in  standard  SPH  stress  calculations  [4]. 


In  additon  and  more  importantly,  however,  the  UC  approach  should  also  help  to 
address  the  tension  instability  present  in  the  conventional  (/?  =  0.5)  SPH.  The  rea¬ 
son  for  this  is  as  follows.  Assume  a  uniform  ID  grid,  as  shown  in  Fig.  1,  with  the 
smoothing  length  2h  =  2L,  where  L  is  the  length  of  each  SPH  element.  For  the 
conventional  approach,  all  nodes  for  the  SPH  gradient  calculations  are  at  least  h  (=L) 
distance  away,  and  node  i  is  not  included  in  any  of  the  calculations  for  stress  and  linear 
momentum  at  i.  In  [2],  it  has  been  shown  that  instability  in  tension  is  governed  by 
(29),  and  for  the  standard  cubic-b-spline  kernel  this  is  satisfied  in  tension  whenever 
x^j  >  0.6h  (approximately).  Unfortunately  for  conventional  SPH  methods,  (29)  is 
satisfied  for  all  the  SPH  nodes  included  in  the  gradient  calculations  at  node  i.  As  men¬ 
tioned  above  for  the  UC  approach,  node  i  is  included  in  the  stress  calculations  at  the 
stress  points  and  (R  <  0.5)  as  shown  in  Fig.  1.  Also,  the  stress  points  and 
ifj  will  be  included  in  the  linear  momentum  calculated  at  node  i.  Thus,  (29)  will  not 
be  satisfied  at  these  two  stress  points.  In  addition,  the  two  SPH  elements  adjacent  to 
i,i  —  I  and  /-Hi,  may  also  contain  stress  points  (depending  on  the  value  of  R)  which 
do  not  satisfy  (29).  Therefore,  the  UC  method  should  have  a  strong  stabilizing  effect 
on  the  SPH  mesh  in  tension,  not  allowing  it  to  prematurely  fragment. 
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6.  NUMERICAL  IMPLEMENTATION  -  SPHID 
For  the  unconventional  (UC)  fonnulation  in  our  ID  code  SPHID,  four  options 
have  been  programmed  for  the  calculation  of  stress  and  linear  momentum.  These  op¬ 
tions  are  specified  by  the  parameter  JSP.  Gradient  calculations  used  to  determine  the 
stress  j  2 )  employ  either  (6)  or  (7)  in  this  work.  Linear  momentum  calculations  are 
based  on  either  (9)  or  (10).  Recall  that  (9)  does  not  conserve  momentum,  but  for  the 
UC  formulation  (6)  and  (9)  are  more  accurate  than  (7)  and  (10).  Equations  (7)  and  (10) 
represent  the  standard  SPH  approach  in  which  additional  approximations  have  been  in¬ 
troduced  [4,5]. 


JSP  =  1  corresponds  to  the  use  of  (7)  and  (10)  or  the  standard  SPH  particle  equa¬ 
tions,  but  with  the  UC  formulation  in  which  the  stress  points  are  not  located  at  the  SPH 
nodes  (0  <  /?  <  0.5).  For  JSP  =  2,  (6)  is  employed  for  the  stresses  while  (9)  is  used 
for  linear  momentum  calculations.  JSP  =  3  uses  (6)  for  the  stresses  and  (10)  for  linear 
momentum.  Thus,  options  JSP  =1  and  3  introduce  additional  approximations  for  the 
UC,  but  (10)  does  enforce  the  conservation  of  momentum.  Also  it  is  noted  that  for  the 
linear  momentum  calculations  using  either  (9)  (JSP  =2)  or  (10)  (JSP  =  1  or  3),  each  of 
the  two  stress  points  within  a  particular  SPH  element  (particle)  is  assumed  to  possess 
one  half  of  the  total  mass  of  that  element  -  this  allows  retention  of  the  particle  concept. 
In  addition,  when  using  (10)  (JSP  =  1  or  3),  the  stress  at  the  typical  SPH  node  /,  since 
it  is  not  calculated,  is  assumed  to  simply  be  the  average  of  the  two  stress  points  and 
ijj  (see  Fig.  1)  in  that  element. 


The  fourth  option  (JSP  =  4)  in  SPHID  uses  the  same  equations  as  JSP  =2  ,  or  (6) 
for  stress  calculations  and  (9)  for  linear  momentum.  However,  with  JSP  =  4  the  stress 
(7  j  j  is  calculated  through  (6)  and  (27)  and  (28)  only  at  the  element  centroid  (as  if 
R  =  0.5  in  Figure  1)  and  it  is  assumed  that: 


U 


'11 


’ll 
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(31) 
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So  G  j  j  is  assumed  to  be  constant  in  SPH  element  i.  Linear  momentum  is  then  calcu¬ 
lated  using  (9)  but  with  the  stresses  applied  not  at  i  (the  centroid)  but  rather  at  and 
as  specified  in  the  UC  approach  when  0<R<  0.5 .  The  option  JSP  =  4,  thus  rep¬ 
resents  a  compromise  between  the  conventional  and  unconventional  SPH  approaches, 
that  may  at  least  help  to  stabilize  the  mesh  in  tension,  although  not  addressing  stress 
oscillation  or  “sawtooth”  behavior. 

Ghost  particles  are  employed  for  the  application  of  essential  boundary  conditions. 
See  [5]  for  the  enforcement  of  free  and  fixed  end  conditions.  Also,  explicit  time  inte¬ 
gration  in  the  form  of  a  modified  central  difference  (MCD)  method  is  applied  to  the 
linear  momentum  equations.  Taylor  and  Flanagan  [9]  briefly  discuss  the  MCD  method 
which  consists  of  a  forward  difference  to  compute  the  velocities  followed  by  a  back¬ 
ward  difference  to  calculate  the  displacements.  In  equation  form  the  MCD  is: 

y(t)  =  y(t-dt)  +dt‘A(t-dt)  (32a) 

U(t)  =  U(t-dt)  +dt-y(t)  (32b) 

where  A  is  the  acceleration  and  dt  is  the  time  increment. 

Overall,  the  SPH  ID  program  is  a  relatively  simple  ID  code,  but  it  will  be  a  useful 
platform  to  demonstrate  our  unconventional  approach.  Finally,  in  SPH  ID  only  a  linear 
elastic  material  model  has  been  implemented  at  this  point. 

7.  APPLICATION 

In  this  section  as  test  of  our  unconventional  approach,  the  SPH  ID  code  will  be 
applied  to  the  elastic  ID  bar  described  in  Fig.  2.  The  bar  is  fixed  at  the  right  end  B  and 
the  left  one  quarter  of  the  bar  is  given  an  initial  velocity  of  Vo=-5  m/sec,  thus  putting 
the  bar  in  tension  initially.  Standard  SPH  methods  (/?  =  0.5)  can  not  solve  this  prob¬ 
lem  due  to  the  tension  instability  that  will  immediately  develop.  As  indicated  in  Fig. 
2a,  the  SPH  grid  is  very  coarse  with  only  40  uniform  SPH  elements  (particles)  used  in 
the  model.  A  comparable  finite  element  model  using  the  ABAQUS  [10,1 1]  program 
that  consists  of  40,  2D  solid  elements  is  described  in  Fig.  2b. 
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Fig.  3  presents  the  displacement  time  history  of  the  left  end  A  (SPH  node  1  actu¬ 
ally)  for  the  UC  results  with  R  =  0.25  (so  the  stress  points  are  at  the  quarter  points 
of  the  SPH  elements).  Also  included  are  the  finite  element  method  (FEM)  results  using 
ABAQUS  with  implicit  time  integration,  which  is  unconditionally  stable.  The  time 
step  used  for  the  explicit  SPH  calculations  as  well  as  for  the  implicit  FEM  model  is  dt 
=  0.4  E-6  sec.  This  time  step  is  based  upon  an  estimate  of  the  Courant  stability  limit 
of  approximately  0.66  E-6  sec.,  which  is  determined  by  dividing  the  length  of  the  typ¬ 
ical  SPH  element  by  the  wave  speed  C. 

As  indicated  in  Fig.  3,  the  SPHID  results  (solid  line)  with  R  =  0.25  are  very 
close  to  the  ABAQUS  predictions  (dashed  line)  for  the  displacement  history  at  the  left 
end  of  the  bar,  point  A.  Some  slight  phase  difference  between  the  two  analyses  devel¬ 
ops  later  but  this  is  probably  due  to  the  use  of  an  explicit  solver  for  SPHID  and  an  im¬ 
plicit  solver  for  ABAQUS,  as  well  as  basic  differences  between  the  SPH  and  FEM 
forms  of  discretization. 

The  JSP  =  2  option  was  used  to  generate  Fig.  3.  This  option  in  SPHID  employs 
(6)  for  stress  calculations  and  (9)  for  linear  momentum.  The  use  of  the  other  options, 
JSP  =1,3  which  employ  (10)  and  do  conserve  momentum  as  well  as  JSP  =  4,  pro¬ 
duced  results  that  eventually  went  unstable.  JSP  =  4  went  unstable  the  quickest  in  the 
analysis.  Reducing  the  timestep  dt  and  increasing  the  artifical  viscosity  did  not  stabi¬ 
lize  the  SPHID  calculations  for  JSP  =1,  3  or  4.  Only  JSP  =2  remained  stable  and  all 
the  results  in  this  section  used  that  option. 

Various  values  for  R  in  the  approximate  range  of  0<R<  0.40  produced  SPH 
results  very  similar  to  Fig.  3  in  which  R  =  0.25 .  For  the  range  0.4  <  /?  <  0.5 ,  the 
response  of  the  bar  tended  to  be  become  unstable  as  the  two  stress  points  approached 
the  centroid  of  the  SPH  element  (i?  =  0.5). 

Fig.  4  indicates  the  predicted  time  history  for  the  velocity  of  the  left  end  A  for 
SPHID  and  ABAQUS.  In  general,  the  agreement  is  excellent  with  again  some  slight 
phase  differences  developing  as  the  analysis  goes  on. 
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In  Figs.  5  and  6,  the  stress  j  j  at  SPH  node  1 1  (stress  point  1 1^)  is  compared  to 
(J  j  j  at  the  centroid  of  finite  element  1 1  -  point  C  in  Figs.  2a  and  2b.  Fig.  5  indicates 
the  extended  time  history  of  O' j  j .  Fig.  6  plots  the  results  only  up  to  .0{X)2  sec.  so  that 
the  differences  in  the  two  analyses  are  more  evident.  As  indicated  in  Fig.  6,  early  in 
the  analysis  especially  the  SPHID  stress  results  do  tend  to  fluctuate  more  than  the 
ABAQUS  results  .  In  general,  the  SPHID  predictions  for  Ojj  compare  reasonably 
well  to  ABAQUS,  but  not  as  well  as  did  the  displacement  and  velocity  at  A  as  indicated 
in  Figures  3  and  4.  This  is  not  too  surprising  since  stress  is  determined  by  differenti¬ 
ating  the  displacement,  and  the  mesh  is  very  coarse  with  only  40  SPH  and  40  FEM 
elements  used  in  the  analyses. 

The  ID  bar  is  next  given  an  initial  velocity  of  Vq  =  5  m/sec,  thus  putting  the  bar 
in  compression  initally.  This  is  done  to  determine  the  effect  of  the  UC  formulation  on 
the  stress  oscillation  or  “sawtooth”  behavior  encountered  in  the  standard  or  conven¬ 
tional  SPH  approach  (R  =  0.5).  Putting  the  bar  initially  in  compression  allows  a 
comparison  between  the  UC  and  the  conventional  SPH,  since  the  latter  is  not  stable  in 
tension.  Fig.  7  shows  a  snapshot  of  the  stress  Gj  j  along  the  entire  bar  at  t  =  1.0  E-5 
sec.  The  UC  results  with  R  =  0.25  are  indicated  by  the  solid  line,  and  the  conven¬ 
tional  SPH  results  in  which  R  =  0.5  are  depicted  by  the  dashed  line.  In  general,  we 
see  that  the  UC  results  are  much  smoother  indicating  a  lessening  of  the  “sawtooth”  be¬ 
havior. 


8.  CONCLUSIONS 

In  this  work,  tension  instability  has  been  addressed  for  SPH  methods.  A  new  un¬ 
conventional  (UC)  approach  has  been  presented  in  which  the  stresses  are  computed  at 
points  away  from  the  SPH  nodes.  The  location  of  these  two  stress  points  within  a  ID 
SPH  element  (particle)  is  controlled  by  the  parameter  R.  The  UC  approach  removed 
the  tension  instability  in  the  bar  considered  for  R  in  the  approximate  range  of 
0<R<  0.45 .  No  optimum  value  of  R  was  found  in  general.  UC  displacement  and 
velocity  results  were  in  excellent  agreement  with  a  comparable  ABAQUS  finite  ele- 
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ment  model.  Axial  stress  comparison  to  the  ABAQUS  model  were  not  quite  as  good, 
but  that  probably  is  to  be  expected  given  the  coarseness  of  the  meshes  and  the  non-lo¬ 
cal  nature  of  the  SPH  method,  which  interpolates  outside  the  SPH  element  (unlike 
FEM).  In  a  comparison  to  standard  SPH  (/?  =  0.5)  for  the  bar  in  compression  ini¬ 
tially,  the  UC  seemed  to  alleviate  the  stress  oscillation  or  “sawtooth”  effect. 

The  UC  approach  was  successful  in  removing  tension  instability,  but  only  for  the 
option  JSP  =  2,  in  which  (6)  and  (9)  were  used  for  the  stress  and  momentum  calcula¬ 
tions.  Equations  (6)  and  (9)  are  not  widely  used  in  the  SPH  literature,  and  (9)  does  not 
conserve  momentum.  These  are  important  points  to  note  with  the  UC  approach. 

Unfortunately,  the  JSP  =4  option  (like  JSP  =1  and  3)  failed  to  remain  stable  for 
the  UC  method.  If  this  option  had  been  successful,  modifications  to  existing  SPH 
codes  would  have  been  easier  and  less  expensive  since  only  one  stress  point  (the  cen¬ 
troid)  within  each  particle  would  have  to  be  tracked. 

Overall,  based  on  the  results  from  the  simple  IDbar,  the  use  of  the  unconventional 
approach  for  SPH  calculations  is  very  encouraging.  This  approach  can  be  extended  to 
2D  and  3D  problems,  but  this  will  be  computationally  expensive  since  additional  stress 
points  within  each  SPH  element  have  to  be  tracked  in  the  analysis  besides  just  the  cen¬ 
troid.  Also,  extension  to  2D  and  3D  may  require  a  rethinking  of  the  overall  particle 
concept.  Finally,  perhaps  a  method  analogous  to  hourglass  control  for  reduced  inte¬ 
gration  in  finite  element  techniques  [11]  may  be  possible  with  the  UC  method. 
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o  SPH  Node  ( at  Centroid) 

X  Stress  Point  ( 0  <  R  <  0.5  ) 
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Figure  1.  Typical  ID  SPH  elements  for  the  unconventional  approach.  When 
R  =  0.5  have  the  conventional  SPH  method. 
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Figure  2.  A  simple  bar  given  an  initial  velocity:  (a)  SPHID  grid,  (b)  ABAQUS 
FEM  grid. 
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Figure  5.  Extended  stress  history  at  node  1 1  (point  C). 


Figure  6.  Early  stress  history  at  node  1 1  (point  C) 


Figure  7.  SPHID  UC  and  conventional  stress  predictions  for  the  bar  at  t  =  l.OE-5 
sec  (Vq  =  .5  m/sec). 
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